﻿using MathNet.Symbolics;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.Statistics;
using MathNet.Symbolics;
using Expr = MathNet.Symbolics.SymbolicExpression;
using MathNet.Numerics.Distributions;

namespace NumericalAnalysis.DataInterpolation
{
    public class NewtonDifferenceQuotient
    {

        static Expr polynomial = 0;   //多项式表达式
        List<double> polyCoefficients = new List<double>(); //各项系数 幂次高至低
        List<int> coefficienOrders = new List<int>(); //各项幂次

        static Vector<double> theVec = Vector<double>.Build.Dense(2);
        List<Vector<double>> vectors = new List<Vector<double>>();


        //public Expr[,] diffQuotMat { get; set; } = null;
        public Matrix<double> diffQuotMat { get; set; }   

        List<double> xs = new List<double>();
        List<double> ys = new List<double>();

        int count = 0;

        /// <summary>
        /// 初始化
        /// </summary>
        /// <param name="_vectors"></param>
        public NewtonDifferenceQuotient(IEnumerable<Vector<double>> _vectors)
        {
            count = _vectors.Count();
            vectors = _vectors.ToList();

            xs = _vectors.Select(v => v[0]).ToList();
            ys = _vectors.Select(v => v[1]).ToList();
        }

        public NewtonDifferenceQuotient(IEnumerable<Tuple<double, double>> _tuples)
        {
            count = _tuples.Count();
            vectors = _tuples.Select(v =>
            {
                var vec = Vector<double>.Build.Dense(2);
                vec[0] = v.Item1;
                vec[1] = v.Item2;
                return vec;
            }).ToList();

            xs = _tuples.Select(v => v.Item1).ToList();
            ys = _tuples.Select(v => v.Item2).ToList();
        }

        public Matrix<double> DiffQuotient ()
        {
            diffQuotMat = Matrix<double>.Build.Dense(count,count,0.0d);

            for (int i = 0; i < count; i++)
            {
                diffQuotMat[i, 0] = ys[i];
            }

            for (int j = 1; j < count; j++)
            {
                for (int i = j; i < count; i++)
                {
                    diffQuotMat[i, j] =
                        (diffQuotMat[i, j - 1] - diffQuotMat[j - 1, j - 1])
                        /
                        (xs[i] - xs[j - 1]);

                    //diffQuotMat[i, j] =
                    //   (diffQuotMat[i, j - 1] - diffQuotMat[i-1, j - 1])
                    //   /
                    //   (vectors[i][0] - vectors[i-j][0]);
                }
            }

            return diffQuotMat;
        }
          
        /// <summary>
        /// 常规公式解法
        /// </summary>
        /// <returns></returns>
        public Expr FitInterpolation()
        {
            var diffQuot = DiffQuotient();
            List<double> dq = diffQuot.Diagonal().ToList();

            var xs = vectors.Select(v => v[0]).ToList();

            //polynomial = 0;
            if (polynomial.ToString().Contains("t")) return polynomial;

            var t = Expr.Variable("t");

            polynomial = dq[0];
            Expr termPoly = t - xs[0] ;
            for (int i = 1; i < count; i++)
            {
                polynomial += dq[i]*termPoly;
                termPoly *= (t - xs[i]);
            }

            //显示多项式
            Console.WriteLine("--------------------------------");
            Console.WriteLine("多项式 polynomial:");
            Console.WriteLine(polynomial);

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式简化 polynomial-ExponentialSimplify:");
            //Console.WriteLine(polynomial.ExponentialSimplify());

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式展开式 polynomial-Expand:");
            //Console.WriteLine(polynomial.Expand());

            polynomial = polynomial.Expand();

            return polynomial;
        }


        /// <summary>
        /// 从给出的向量计算
        /// </summary>
        /// <param name="sourceVecs"></param>
        /// <returns></returns>
        public IList<Vector<double>> CalculateBy(IList<Vector<double>> sourceVecs)
        {
            for (int i = 0; i < sourceVecs.Count(); i++)
            {
                double xValue = sourceVecs[i][0];
                var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
                double yValue = polynomial.Evaluate(symbols).RealValue;
                //double   yValue = CalculateBy(xValue);
                sourceVecs[i][1] = yValue;
            }
            return sourceVecs;
        }

        /// <summary>
        /// 从给出的数列计算
        /// </summary>
        /// <param name="xValues"></param>
        /// <returns></returns>
        public IList<double> CalculateBy(IList<double> xValues)
        {
            List<double> yValues = new List<double>();
            for (int i = 0; i < xValues.Count(); i++)
            {
                var symbols = new Dictionary<string, FloatingPoint> { { "t", xValues[i] } };
                double yValue = polynomial.Evaluate(symbols).RealValue;

                yValues.Add(yValue);
            }
            return yValues;
        }

        /// <summary>
        /// 从给出的数计算
        /// </summary>
        /// <param name="xValue"></param>
        /// <returns></returns>
        private double CalculateBy(double xValue)
        {
            //复杂多项式避免使用Compile方法计算  将严重拖慢性能
            //Func<double, double> func = polynomial.Compile("t"); 
            //double yValue = func(xValue);

            var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
            double yValue = polynomial.Evaluate(symbols).RealValue;

            return yValue;
        }



    }



    public class NewtonDifferenceQuotientTest
    {


        /// <summary>
        /// 测试
        /// </summary>
        /// <param name="exsistPnts">存在插值点</param>
        /// <param name="needPntXs">需要求的点自变量</param>
        /// <returns></returns>
        public static IList<Tuple<double, double>> NewtonDifferenceQuotientCalc(IList<Tuple<double, double>> exsistPnts, IList<double> needPntXs)
        {
            var t = Expr.Variable("t");
            NewtonDifferenceQuotient interp = new NewtonDifferenceQuotient(exsistPnts);
            var expression = interp.FitInterpolation(); 

            var needPntYs = interp.CalculateBy(needPntXs);

            var tuples = new List<Tuple<double, double>>();
            for (int i = 0; i < needPntXs.Count; i++)
            {
                var tuple = new Tuple<double, double>(needPntXs[i], needPntYs[i]);
                tuples.Add(tuple);
            }
            return tuples;
        }


        public static (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleHomeWork()
        {
            var t = Expr.Variable("t");
            List<Tuple<double, double>> extVecs = new List<Tuple<double, double>>();

            List<double> values = new List<double>()
                { 12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 };
            for (int i = 0; i < values.Count; i++)
            {
                double xk = i * 2;
                double yk = values[i];
                Tuple<double, double> tuple = new Tuple<double, double>(xk, yk);
                extVecs.Add(tuple);
            }
            var interpYs = NewtonDifferenceQuotientCalc(extVecs, extVecs.Select(v => v.Item1).ToList()).ToList();

            List<double> array2 = NumericalService.LinSpace(0 + 0.1, 24 - 0.1, 0.1);
            var checkdYs = NewtonDifferenceQuotientCalc(extVecs, array2).ToList();

            return (interpYs, checkdYs);
        }


        /// <summary>
        /// sin函数
        /// </summary>
        /// <returns></returns>
        public static (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleSin()
        {
            var t = Expr.Variable("t");
            List<Tuple<double, double>> extVecs = new List<Tuple<double, double>>();
            //构造sin函数插值点
            List<double> xs = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.1);
            extVecs = NumericalService.CalcByFunc((x) => { return Math.Sin(x); }, xs).ToList();

            ////进行插值 并 返回检查点    这里插值点就是检查点  
            var interpYs = NewtonDifferenceQuotientCalc(extVecs, extVecs.Select(v => v.Item1).ToList()).ToList();

            List<double> array2 = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.02);
            var checkdYs = NewtonDifferenceQuotientCalc(extVecs, array2).ToList();

            return (interpYs, checkdYs);
        }


        /// <summary>
        /// 龙格函数
        /// </summary>
        /// <returns></returns>
        public static (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleRunge()
        {
            var t = Expr.Variable("t");
            List<Tuple<double, double>> extVecs = new List<Tuple<double, double>>();
            //构造runge函数插值点
            Func<double, double> rungeFunc = (x) => { return 1 / (x * x + 1); };
            List<double> array = NumericalService.LinSpace(-5, 5, 0.5);
            extVecs = NumericalService.CalcByFunc(rungeFunc, array).ToList();
            ////进行插值 
            var interpYs = NewtonDifferenceQuotientCalc(extVecs, extVecs.Select(v => v.Item1).ToList()).ToList();

            List<double> array2 = NumericalService.LinSpace(-4, 4, 0.02);
            var checkdYs = NewtonDifferenceQuotientCalc(extVecs, array2).ToList();

            return (interpYs, checkdYs);
        }

    }

}
